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Abstract 

Numerical data structures for positive dimensional solution sets of polynomial systems are 
sets of generic points cut out by random planes of complimentary dimension. We may represent 
the linear spaces defined by those planes either by explicit linear equations or in parametric 
form. These descriptions are respectively called extrinsic and intrinsic representations. While 
intrinsic representations lower the cost of the linear algebra operations, we observe worse 
condition numbers. In this paper we describe the local adaptation of intrinsic coordinates 
to improve the numerical conditioning of sampling algebraic sets. Local intrinsic coordinates 
also lead to a better stepsize control. We illustrate our results with Maple experiments and 
computations with PHCpack on some benchmark polynomial systems. 
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1 Motivation, Definitions, and Problem Statement 

A polynomial system /(x) = 0, x = (xi,X2, . . . ,Xn), defines an algebraic set /"^(O) C C". The 
polynomials of / belong to C[x]. We assume (for simplicity of exposition throughout the paper): 

1. /^^(O) is pure dimensional, k is its codimension, so dim(/^-'^(0)) = n — k; 

2. /(x) = is a complete intersection, and in particular: / = (/i, /2, . . . , fk); 

3. /^^(O) is reduced, i.e.: of multiplicity one. 

To remove the third assumption, a deflation operator [21] (see also [10]) as proposed in 
§13.3.2] should be applied. The first two assumptions are made for notational convenience. 

The numerical treatment of positive dimensional algebraic sets was first proposed in 
and elaborated in a series of papers by the authors of [32] and the second author, see also [30l 
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for another introduction. The algorithms in numerical algebraic geometry are implemented in 
PHCpack [33] and Bertini [3] (see [6] and [26]) and can be executed via MATLAB (or Octave) [15] . 
Maple [20], and Macaulay 2 [l9]. 

One of our benchmark examples is a family of systems, defined by all adjacent minors of a 
general 2-by-3 matrix ([12], [IB]): 



/(x) 



Xii Xi2 Xi3 
X21 X22 X23 

For this example, we have n = 6, k - 
n — k = 4. To compute deg(/~^(0)), we add n 



X\\X22 - X2\X\2 = 
X12X22. - X22Xvi = 0. 

2, and we have a complete intersection: dim(/~^(0)) 
k general linear equations L(x) = to /(x) = 



(1) 







and solve {/(x) = 0, L(x) = 0}. Generic points on the solution set defined by the system for all 
adjacent minors of a general 2-by-3 matrix satisfy (for random coefficients Qj E C): 



(2) 



X11X22 - X21X12 = 
X12X22. - X22X13 = 
cio + ciixii + C12X12 + C13X13 + C14X21 + C15X22 + C16X23 = 

C20 + C21X11 + C22X12 + C23a:i3 + C24X2I + C2r>X22 + C26X23 = 
C30 + C31X11 + C32X12 + C33X13 + C34X21 + C35X22 + 036X23 = 
, C40 + C41X11 + C42X12 + C43X13 + C44X21 + C45X22 + C46X23 = 0. 

Except for an algebraic set in the coefficient space Cij for L, the system above has four solutions, 
we have four generic points for all adjacent minors of a general 2-by-3 matrix, so deg(/~^(0)) = 4. 

To save work, reducing the number of variables from 6 to 2, we choose a different representation 
for the linear space defined by the equations iv(x) = 0, representing the 2-plane L~^(0) in as 



(3) 



spanned by an offset point b G and an orthonormal basis {vi, V2}. The tuple (^1,1^2) defines 
intrinsic coordinates for the generic points, introduced in [29j to speedup the algorithms of [28] . 

The reduction from six to two variables reduces the cost of solving linear systems by a factor of 
nine. This reduction improves the efficiency of Newton's method when computing sample points 
on the algebraic set, one of the basic operations in numerical algebraic geometry [32] • 

For any f^^{0) G C" with dim(/^^(0)) = n — k, we use a general A;-plane L to compute 
generic points. This general A;-plane L may be defined in two equivalent ways: 



" xu ' 




' h ' 




" Vll ' 




" V21 ' 
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1. L{x) = is a system of n — k general linear equations in x, 

Vfcl G C"^^ with V*V 



2. b G C" is an offset point, and V 
orthonormaQ basis of vectors. 



Vl V2 



Ik, i.e.: V is an 



^Although it suffices to require that the columns of the matrix V are hnearly independent, the orthonormahty 
condition V*V — Ik (using complex conjugated inner products and Ik is the fc-by-fc identity matrix) is beneficial. 
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If linear equations define L, solving {/(x) = 0, L(x) =0} gives generic points in their usual form 
that we call an extrinsic coordinate representation. Using (b, V) for L gives intrinsic coordinates 
^ = (6, 6, • • • , 6) for generic points x: 



X = b + ClVi + + • • • + CfcVfc = h + Vt 



(4) 



With intrinsic coordinates for generic points, the original variables x become place holders when 
solving /(x = b + V^) = 0. In Figure [H we outline the two ways to compute generic points. 



L- 



Ke 



Ki 



X 



Figure 1: A commutative diagram for extrinsic x and intrinsic ^ coordinates of generic points. 
The vertical arrows require linear algebra while the horizontal arrows involve the solution of 
polynomial systems. Their sensitivities are determined by condition numbers Ke and K^. 



In shorthand notation, the general /c-plane L is represented as (b, V) and we use intrinsic 
coordinates ^ € to denote the generic points. When sampling points, the moving L from 
(b, V) to (c, W), is done via the obvious homotopy: 

/x= (l-t)b + tc + [{l-t)V + tW) o_ .5) 

y moving offset point moving basis vectors J ' 

As t moves from to 1, the solution paths ^(t) are tracked with predictor-corrector methods and 
give new generic points on /~^(0). For introductions to path following and continuation methods 
we refer to [1] and [23|, see also [22] . 

While the diagram in Figure [T] commutes for exact operations, using floating-point arithmetic 
forces us to take into account condition numbers. These condition numbers bound the growth of 
the relative errors on the solutions as a consequence of relative errors on the input data. As we 
keep / fixed during the computation we only consider relative errors on the representations of the 
/c-plane L. Formally, we introduce condition numbers Ke and Kj on the extrinsic and intrinsic 
coordinate representations respectively as 

[lAxll IIALII , IIA^II ||A(b,n|| 

Going to intrinsic coordinates, we observe a worsening of the numerical conditioning: Kj ^ Ke- 
Note that the original problem is well conditioned, in other words, Ke is expected to remain 
small, because random choices for L avoids places where the Jacobian matrix of / drops rank. 

To get a first intuition why Kj ^ Ke, consider the binomial expansion of a monomial x^x^ 
of /. If we evaluate x^.^x^' at xi = 61 + ^\V\ and X2 = 62 + ^2^2-, we compute 



(61 + iiv^r ib2 + ^2V2r = E 7 E ; 4i^2V2r-^ (7) 
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and we see that any sparse structure of / will be destroyed. Moreover, the binomial coefficients 
in ([7]) inflate the variation among the coefficients in /. 

In general, we may write /(b + + A^)) = /(b + V^) + A/. The trouble is that, even for 
small ||A^|| we may experience very large ||A/||. 

While using multiprecision arithmetic during path tracking [5] may avoid these numerical 
instabilities, using multiprecision numbers significantly slows down the computations and when 
the coefficients are known with limited accuracy, applying multiprecision arithmetic may give 
misleading answers. Better stepsize control strategies [7j will also be effective for our problems, but 
like in dealing with the high powers of the continuation parameter of polyhedral homotopies |18j . 
our approach in this paper is specific to the type of homotopies. To deal with the numerical 
instabilities of using intrinsic coordinates, we propose the use of local intrinsic coordinates. We 
define local coordinates in the next section. In section 3, we present an algorithm to track 
a solution path using intrinsic coordinates, along with an a priori stepsize control evaluation 
strategy. Computational results are discussed in the section 4. 

Acknowledgement. We thank Professor Hiroshi Murakami for his remarks made after the 
presentation of the first author at the session of Symbolic and Numeric Computation at ACA 
2009. His remarks led us to local intrinsic coordinates. 

2 Local Intrinsic Coordinates 

In this section we define local intrinsic coordinate representations of generic points and address 
the improved numerical conditioning. 

What if we could keep ||^|| small? Writing Greek symbols badly, the ^ looks close enough to 
an epsilon, and then reconsidering the binomial expansions in ([7]): 

ih + ^ivir {b2 + i2V2r = {bT + aibl^-'CiVi + 0(e?)) + a2b-,'-'^2V2 + Ofe')) (8) 

If we assume that is infinitesimally small, then we ignore the second order terms 0{(,f ,(,1^2,^,2) ■ 
For general polynomials /, writing /(b + V$,) = /(b) + A/, the omission of the higher order 
terms leads to: ||A/|| is 0(||y^||). Because we may select for the orthonormal basis V a nice 
numerical representation, we have that 0(||y^||) is 0(||^||) and therefore: ||A/|| is ||0(^)||. 

To keep ||^|| small, we now propose to use the extrinsic coordinates of the generic point as the 
offset point for a fc-plane. In particular, for d = deg(/~^(0)) and d generic points {zi, zi, . . . , z^} 
on f^^{0), consider: 

x = z^ + l^|, £=1,2, (10) 

Because L(z^) = for all generic points, all zi + represent the same /c-plane L. Given an 
orthonormal basis V for a /c-plane and a set {zi,zi, . . . ,Zd} of d generic points on /~^(0), the 
local intrinsic coordinates to represent /~^(0) are defined by the tuple ({zi, zi, . . . , z^}, V). 

Obviously, the transition from global intrinsic coordinates x = b + to local intrinsic coor- 
dinates is performed by a mere evaluation of b + V^. The close relation between local intrinsic 
coordinates and extrinsic coordinates will yield improved condition numbers. 
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To define the condition number Ke of a zero z of F := {/(x) = 0,L(x) = 0}, we consider 
the application of Newton's method: 

F'(z)Az = -F(z), Ke:=k{A), (11) 
=A 

where F' is the matrix of all partial derivatives of F and n{A) is the condition number of the 
Jacobian matrix A of F at z. Because we assume that /(x) = is a complete intersection, 
AAz = —F{z) is a well defined n-by-n linear system. Strictly speaking, as we keep / fixed and 
vary only the linear equations L(x) = 0, we will have Ke < but because generic points are 

always well conditioned this distinction is very minor. 

In local intrinsic coordinates we replace x by z + and the application of Newton's method 
leads to 

f'{z + Vi)A^ = -fiz + VO, Ku:=k{B), (12) 

=B 

where /' is the matrix of all partial derivatives of / and k{B) is the condition number of the 
Jacobian matrix of / at ^. Because we assume that /(x) = is a complete intersection, 
B/S.^ = — /(^) is a well defined k-hy-k linear system. We define k{B) as Klj, the condition 
number of z represented in local intrinsic coordinates. 

Observe the similarity of (jlip with (jl2p as we write (jlip more explicitly as 



L'fz) 



Az 



/(z) 
L(z) 



(13) 



where L' contains all partial derivatives of L. 

Proposition 2.1 With Ke and K^j as defined in (jlip and (jl2p respectively: K^j ~ Ke- 

Proof. Because in local intrinsic coordinates: ^ = 0, it does no longer make sense to consider 
relative errors. Moreover, without loss of generality we may always choose coefficients of the 
planes so that ||L|| = 1 and \ \V\\ = 1. Using homogeneous coordinates for /, we assume we work 
in an appropriate affine chart so that also ||z|| = 1. Then the meaning for the condition numbers 
Ke and Kli are in the inequalities 

[|Azl| < ifij||AL|| and [|A^|| < KL/l|A(b, F)||, (14) 

where b is an offset point and V are directions in a parametric representation of L. 
Using the commutative diagram of Figure [TJ we relate Az and A^: 

z + U(^ + A^) = z + Az ^ Az = UA^ as ^ = 0. (15) 

Looking at norms: 

||Az|| = IIUA^II = IIA^II as||U|| = l. (16) 

Because V*V = Ik, all eigenvalues of V lie on the complex unit circle and multiplication with V 
is norm preserving. 

In local intrinsic coordinates, for ^ = 0, changes AV in the orientation of L do not influence ^. 
So we have A(b,U) = Ab and in case we may consider ||AL|[ ~ ||A(b,U)|[. Thus in (jl4p we 
may interchange Ke with K^j, so K^j Ke- □ 
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3 A Rescaling Algorithm 



In this section we consider the sampHng of algebraic sets using local intrinsic coordinates. We 
define a rescaling algorithm and address its numerical stability. In addition, using local intrinsic 
coordinates leads to a better stepsize control. 

Generic points {zi, zi, . . . , z^} are offset points for a A:-plane L with directions in the orthonor- 
mal matrix V. In local intrinsic coordinates, moving from (z£,V) to (b,VF), as t goes from 
to 1, the deformations are defined by 

/(x = (i-t)z, + tb + iyo = 0. (17) 

In contrast to the obvious homotopy in ([5]), we see that only the offset point moves. We immedi- 
ately switched from the current directions in V to the new orthonormal basis W because ^ = 
in local intrinsic coordinates. But this is only a first indication of the potential of working with 
local intrinsic coordinates, we can do better than ()17p . 

Instead of using (jl7p and moving to b, we point out that any point in the fc-plane L can serve 
as an offset point. Therefore, we should choose the best offset point, i.e.: the point closest to the 
current generic point. To compute the closest point, let c be the orthogonal projection of z^ onto 
the A:-plane L. For some step size h, we then consider: 

/ (x = z, + hie - ze) + W0 = (18) 

and apply Newton's method to find the correction A^, as illustrated in Figured 




Figure 2: Schematic of one predictor-corrector step of the new sampling algorithm, moving from 
the point z to the point where L meets /^^(O). The line L is defined by b + ,^w. Using step size 
h, the prediction h{c — z) added to z occurs in a direction orthogonal to L while the correction 
A^w is parallel to L. 

After each step, we add the correction term (A^w in Figure[2|) to the offset point, rescaling the 
intrinsic coordinates to local intrinsic coordinates at the end of the correction stage. Pseudocode 
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for one predictor-corrector step is given in Algorithm I3.H going from one generic point z E 
/~^(0) n K, where K is the current fc-plane, towards L the target A;-plane. 

Algorithm 3.1 (one predictor-corrector step in local intrinsic coordinates) 

Input: / = (/i, /2, . . . , A), /i(x) G C[x], z = 1, 2, . . . , fe; dim(/-i(0)) =n-k 

b G C"; offset point of k-plane L 

W = [wi W2 • • • Wjt] G C"^'^, W*W = Ik] orthonormal basis for L 

z G C": /(z) = 0, K{z) = 0; generic point on k-plane K 

h > 0; step size 

e > 0. accuracy requirement 

Output: z, /(z) = 0, L{z) =0: ||z — b|| < ||z — b||. generic point closer to L 



1. 


V := z — b; 

k 


go towards offset point 


2. 


V := V - '^{Wfv)wi; 


move perpendicular to L 


3. 


i=l 

I '~ 


normalize so v = 1 


4. 


z := z + /i v; 


prediction for new generic point 


5. 


z := z; ^ := 0; 


initialize for Newton corrector 


6. 


while 11/(2 + T^^) II > e do 


as long as not accurate enough 


6.1 


■.= f{z + WO/f{z + W^y, 


solve a linear system for 


6.2 


^■.= ^ + At 


update correction 


7. 


z ■.= z + W$. 


rescale to local coordinates 



The orthonormality condition W*W = Ik is important for instruction 2 in the Algorithm 13.11 
because we can compute the projection just via inner products. The number of arithmetical 
operations needed to carry out instruction 2 in Algorithm 13.11 is 0{kn). Without the condi- 
tion W*W = Ifc, this cost (e.g. via Gram-Schmidt orthogonalization) would be at least 0{kn^). 

For the numerical stability of Algorithm 13. 11 we first discuss the relationship between the step 
size h and the accuracy requirement e. If on the one hand h is too small, then the condition 
||/(z + 1^^)11 > e in the while-do instruction 6 of Algorithm 13.11 is directly satisfied. On the 
other hand, if h is too large, satisfying the accuracy requirement of instruction 6 may require 
too many iterations, or Newton's method may not converge at all. We point out that the cost of 
instruction 6.1 is 0{k^) and if / is sufficiently sparse (if evaluation and differentiation go fast), 
then the cost of execution of Newton's method dominates the cost of Algorithm 13. li 

In general path tracking algorithms, the step size h is determined via a feedback mechanism. 
If Newton's method does not converge fast enough, then the step size is reduced. If Newton's 
method needs only two steps or less, then the step size might be enlarged. See [7] for stepsize 
control strategies. The problem with this feedback mechanism is that it comes at the great 
expense of the most costly portion of the predictor-corrector method, i.e.: each reduction of h 
comes at the expense of a failed and thus wasted Newton step. With local intrinsic coordinates, 
we can predict the fitness of the step size with a simple evaluation. For some step size h and 
direction v, we evaluate and estimate the residual as 

||/(x = ze + hw)\\ is ||/(z^) + 0{h)\\ is 0{h). (19) 
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For example, if h = 10~^ and we see that the residual is O(10~^), then it is fair to expect that 
after one iteration of Newton's method, the residual becomes 0(10^^), and then O(10~^) after 
the second iteration. 

In Algorithm 13.21 we define how to cut back on the step size just by evaluation, before the 
start of the Newton correction. 



Algorithm 3.2 (a priori stepsize control by evaluation) 

Input: /= (/i,/2,...,/fc),/i(x) eC[x] ^= 1,2,...,A;; dim(/-i(0)) = n - A: 
z G C": /(z) = 0, K{z) = 0; generic point on k-plane K 

V € C", [|v[| = 1; direction vector 

h > 0; current step size 

5 > 0. threshold to reduce h 

1 > p > 0. reduction factor for h 

Output: h > 0. updated step size 

1- y + ^v)||; evaluate the predicted point 

2. if y/h > 5 then h := ph. reduce the step size 

The reduction of the step size in instruction 2 of Algorithm 13.21 could be followed by another 
evaluation of / to see if y is reduced sufficiently or has become even too small. 

Algorithm 13.21 is called after instruction 3 of Algorithm 13.11 

By application of Algorithm 13.21 occurrences of a diverging Newton's method can be greatly 
reduced because the size of the residual ||/(x = z + P1^^)|| is 0(h). 

We conclude with a quick cost estimate for the total number of Newton steps along one path. 
In sampling for generic points, we typically choose the new random coefficients for the /c-plane as 
complex numbers on the unit circle, so the distance between two /c-planes (and in particular their 
offset points) is 0(1). For h: < /i < 1, we can see that the total number of Newton iterations 
along a solution path is proportional to 1/h. For example if h = 0.01 and we need about 2 or 3 
Newton iterations per step, then the total number of Newton iterations along a solution path will 
vary between 200 and 300. 

The homotopy continuation methods of this paper are different from the so-called linear 
homotopies for which an experimental study to certify path tracking recently appeared in [8]. A 
potential future research direction could be to expand the quick cost estimate of the previous 
paragraph into a formal complexity study, along the lines of ^Oj and [25] . 



4 Computational Results 

Local intrinsic coordinates are available in version 2.3.53 of PHCpack [53]. We first describe 
numerical experiments done with Maple to compare condition numbers of generic points on a 
hypersurface of polynomials of increasing degrees. Then we report preliminary results on small 
benchmark problems with the sampling routines in PHCpack. All computations were done on 
one core of a Mac OS X 3.2 Ghz Intel Xeon. 
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4.1 Condition Number Estimates 



In this section we illustrate the worsening of the conditioning of using global intrinsic coordinates 
on one sparse polynomial. We give data on sampling with zero and nonzero offset vectors and 
relate this experiment to using local intrinsic coordinates. 

To estimate the condition numbers we use LinearAlgebra[EigenConditionNumbers] of 
Maple 12, with UseHardwareFloats set to true, see \23\ Chapter 4]. The corresponding doc- 
umentation pages of Maple 12 refer to p]. For an introduction to the perturbation theory of 
eigenvalues, see e.g.: [Hj §4.3]. 

We consider one sparse polynomial / in n = 10 variables, of increasing degrees d, starting 
with t terms. In addition, we add all the linear terms CjXj, i = 1,2, ... ,n, to avoid ending up 
with the origin as a multiple root. The coefficients are taken on the complex unit circle. The 
particular Maple commands used to generate an / are 

[> n := 10: d := 10: t := 5: 

[> c := -> exp(I*stats[random,uniform[0,2*Pi]] (1)) : 

[> X := [seq(x[i] ,i=l. .n)] : 

[> f := X[l]"d + randpolyCX, coef f s=c ,degree=d-l ,terms=5) + sum(c () *x [i] , i=l . .n) ; 

The first term of f ensures that we have a monic polynomial after substitution /(v^), for vi = 1. 
That / is monic is convenient for the connection with the companion matrix when we look at the 
condition numbers of the corresponding eigenvalue problem. 

To introduce the idea of using different coordinate systems, we respectively use 

x = b + v^ and x = v^, b,vGC'', (20) 

where all coefficients in the vectors are also taken on the complex unit circle. With /(v^) = 
we obtain still a sparse polynomial with all coefficients on the complex unit circle, which is not 
the case with /(b + v^) = 0. The offset vector of b + is responsible for the variation in the 
coefficients and the fluctuation of the condition numbers we observe in our numerical experiments, 
summarized in Table [H 



degrees 




/(b + 


vO =0 




/K) 


= 


ratios of 


ratios of 


of/ 




largest 


smallest 




largest 


smallest 


smallest 


largest 


10 


5 


.91e-01 


9.02e-02 


8 


.81e-01 


4.01e-01 


6.55e+00 


2.20e+00 


20 


2 


.77e-01 


1.76e-03 


8 


.91e-01 


3.31e-01 


1.57e+02 


2.70e+00 


30 


2 


.75e-01 


6. 16e-05 


9 


.49e-01 


7.25e-02 


4.47e+03 


1.31e+01 


40 


4 


.53e-01 


7. 14e-06 


9 


.69e-01 


1.87e-01 


6.34e+04 


5. 17e+00 



Table 1: Estimates for the inverse condition numbers of eigenvalues of the companion matrices of 
/(b + v,^) = and /(v^) = 0. For degrees d = 10, 20, 30, and 40, we list the largest and smallest 
inverse condition numbers. For /(v^) = 0, we see the range between smallest and largest not 
widen that much, whereas for /(b + v^) = 0, the conditioning steadily worsens. 

As we see from Table [H all roots of /(v^) = are well conditioned. To compare the con- 
ditioning of local intrinsic coordinates, we take the first root zi of /(v,^) = and consider the 
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companion matrix A of f{zi + v^) = 0. For increasing degrees, the condition number for the zero 
^ = corresponding to zi is always reported as 1 . OOe+00. The smallest inverse condition num- 
bers of the eigenvalues of A for degrees d = 10, 20, and 30 are respectively 8.42e-05, 1.08e-12, 
and 3.69e-14. This implies that for d = 30 we have lost all accuracy as our working precision 
are the standard hardware floats. 

In this simple Maple experiment we illustrate that, although sampling a hypersurface is re- 
duced to solving univariate polynomial equations, for hypersurfaces defined by polynomials of 
high degrees we cannot use the same representation of a general line to define generic points. If 
we adapt the offset point and switch to local intrinsic coordinates, then the generic points are 
well conditioned. 

4.2 Sampling Benchmark Systems 

The input to the sampling problem is one set of generic points on /^^(O)nL and a new fc-plane K. 
On output is a new set of generic points on /~^(0) fl K. 

The polynomial systems we selected occur in the literature. We briefiy summarize the main 
characteristics of these systems: 

1. All adjacent minors of a general 2-by-n matrix, n = 3,4,.... This is a family of nice 
quadratic equations arising in algebraic statistics [12] . 

2. The cyclic n-roots systems are well known academic benchmarks. If n has a quadratic 
divisor, then the system has a positive dimensional solution set [3]. In our experiments we 
use the cyclic 8-roots system, which has a one dimensional solution set of degree 144. 

3. Griffis-Duffy platforms |13j are architecturally singular mechanisms [17j, their motion cor- 
respond to curves of degree 40 in 8-space [27]. 

For the purposes of this paper, the computation of the first set of generic points is considered as 
given, typically in extrinsic coordinate representation. 

Except for the adjacent minors, the systems are not complete intersections. For m > k, to 
make an m-by-fc system / square, we generate a random k-hy-m matrix C and work with C x f. 

To test the improvement from using local intrinsic coordinates, we sample new generic points 
from the solution sets. Our computational experimental setup consists of three stages: 

(1) Given one set of generic points, we generate another random /c-plane L. 

(2) We then move the given set of generic points to lie on L. 

(3) At the end we check results for accuracy, count #predictor-corrector steps, record elapsed 
cpu times. 

Note that the recorded cpu times are only meant to give an indication on the relative practical 
difficulties of these problems. More relevant are the number of iterations performed by Newton's 
method along the paths. 

In Table [2] we summarize the results. Even as the systems we selected as benchmark examples 
are not challenging, we observe a clear benefit of using local intrinsic coordinates, even for the 
systems defined by quadratic equations. The benefit is perhaps most significant for the cyclic 
8-roots problem where the degree of the ith polynomial equals i. 



10 



polynomial system 


n 


n — k 


d 


^iterations 


timings 


Griffis-Duffy platform 


8 


1 


40 


207/164 


550/535 yusec 


cyclic 8-roots system 


8 


1 


144 


319/174 


5.3/3.2 sec 


all adjacent minors 


22 


12 


1,024 


285/219 


44.6/40.3 sec 



Table 2: Preliminary experiments on three systems. For each system we respectively list the 
ambient dimension n, the dimension n — k of the solution set, and the degree d of the set. We list 
the average number of Newton iterations along a path for intrinsic and local intrinsic coordinates, 
as well as user cpu timings. 

5 Conclusions 

We list at least three advantages of using local intrinsic coordinates for sampling: (1) only the 
offset point moves; (2) the sparse structure of the polynomials is kept; and (3) we can control the 
step size by evaluation. Applications to numerical algebraic geometry include (1) implicitization 
via interpolation; (2) monodromy breakup algorithm; and (3) diagonal homotopies. In particular, 
local intrinsic coordinates will add to the robustness of our parallel subsystem-by-subsystem 
solver [14j . 
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